Phase Separation of Charge-Stabilized Colloids: 
A Gibbs Ensemble Monte Carlo Simulation Study 
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Fluid phase behavior of charge-stabilized colloidal suspensions is explored by applying a new vari- 
ant of the Gibbs ensemble Monte Carlo simulation method to a coarse-grained one-component model 
with implicit microions and solvent. The simulations take as input linear-response approximations 
for effective electrostatic interactions - hard-sphere- Yukawa pair potential and one-body volume en- 
ergy. The conventional Gibbs ensemble trial moves are supplemented by exchange of (implicit) salt 
between coexisting phases, with acceptance probabilities influenced by the state dependence of the 
effective interactions. Compared with large-scale simulations of the primitive model, with explicit 
microions, our computationally practical simulations of the one-component model closely match the 
pressures and pair distribution functions at moderate electrostatic couplings. For macroion valences 
and couplings within the linear-response regime, deionized aqueous suspensions with monovalent 
microions exhibit separation into macroion-rich and macroion-poor fluid phases below a critical 
salt concentration. The resulting pressures and phase diagrams are in excellent agreement with 
predictions of a variational free energy theory based on the same model. 

PACS numbers: 05.20.Jj, 82.70.Dd, 82.45.-h 



I. INTRODUCTION 



Charge-stabilized colloidal suspensions [1] con- 
taining monovalent microions reportedly can display 
unusual thermodynamic behavior when strongly 
deionized. Puzzling experimental observations in- 
clude liquid- vapor coexistence [2], stable voids [3- 
6], contracted crystal lattices [6-8], and metastable 
crystallites [9] . Such phenomena reveal an extraordi- 
nary cohesion between like-charged macroions that 
appears inconsistent with the purely repulsive elec- 
trostatic pair interactions predicted by the classic 
theory of Derjaguin, Landau, Verwey, and Overbeek 
(DLVO) [10, 11]. Failure of the DLVO theory to 
account for anomalous phase behavior of deionized 
suspensions has prompted many theoretical and sim- 
ulation studies. 

Predictions of a spinodal instability in deion- 
ized charged colloids follow from classical density- 
functional [12-15], extended Debye-Huckel [16-18], 
and linear-response [19-23] theories of a coarse- 
grained one-component model. The predicted phase 
separation is driven by the state dependence of the 
effective electrostatic interactions, including a one- 
body volume energy [24-26]. Such predictions have 
been challenged on grounds that underlying lin- 
earization approximations may fail to describe non- 
linear microion screening [27-29] and neglect strong 
counterfoil association that may renormalize the ef- 
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fective macroion charge [30-34] . The debate is some- 
what complicated, however, by the proximity of the 
unstable parameter regime to the threshold for sig- 
nificant nonlinearity and charge renormalization. 

Some simulations of the primitive model [35, 36], 
with explicit microions interacting via long-ranged 
Coulomb potentials, exhibit clustering of macroions 
at strong electrostatic couplings. Such computation- 
ally intensive simulations become increasingly de- 
manding, however, upon approaching the size and 
charge asymmetries required to directly test predic- 
tions, even when sophisticated cluster moves are in- 
cluded [37]. Therefore, reconciling theories, simula- 
tions, and experiments to clarify the phase behavior 
of deionized charged colloids calls for novel simula- 
tion methods tailored to mesoscale models. 

The main purpose of the present work is to pro- 
pose a new variant of the Gibbs ensemble Monte 
Carlo method suited to modeling density-dependent 
effective electrostatic interactions. As a demonstra- 
tion, we apply the method to deionized charged col- 
loids to test predictions of phase instability. After 
first defining the model system and one-component 
mapping in Sec. II, we briefly summarize the linear- 
response theory of effective interactions and a vari- 
ational free energy theory in Sec. III. The Monte 
Carlo algorithm is next outlined in Sec. IV. Simu- 
lation results are presented in Sec. V, with diagnos- 
tic details deferred to the Appendix. Comparisons 
with theory and primitive model simulations confirm 
previous predictions and illustrate computational 
advantages and limitations of the one-component 
model. Finally, Sec. VI summarizes our conclusions. 
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II. MODEL 



III. THEORY 



A. Primitive Model 

As underlying microscopic model, we adopt the 
primitive model of charged colloids [38]: macroions 
and microions dispersed in a continuum solvent of 
dielectric constant e in a closed volume V. The 
macroions are modeled as charged hard spheres, 
monodisperse in radius a and effective valence Z 
(charge — Ze), and the microions (counterions and 
salt ions) as point charges of valence z. Here we 
assume monovalent microions (z = 1) dispersed in 
water at temperature T — 293 K, corresponding to 
a Bjerrum length \ B = e 2 /(ek B T) = 0.72 nm. As- 
suming N m macroions and N s pairs of dissociated 
salt ions, we have N+ = (Z/z)N m + N s positive and 
N- = N s negative microions. 



A. Linear-Response Theory 

Statistical mechanical descriptions of effec- 
tive electrostatic interactions, including density- 
functional [12-15], extended Debye-Huckel [16-18], 
and response [19-22] theories, typically invoke lin- 
earization and mean-held approximations for the mi- 
croion free energy [Eq. (3)]. Response theory 
describes the perturbation of the microion densities 
by the "external" macroion electrostatic potential. 
Taking as the unperturbed reference system a uni- 
form gas of microions in the free volume outside the 
macroion hard cores, the microion free energy can 
be expressed as 



F„ = Fr 



plasma 



dA (H rnp ) x — E\ 



(4) 



B. Coarse-Grained One-Component Model 

Long-ranged Coulomb interactions and high 
charge asymmetries render large-scale simulations 
of the primitive model computationally challeng- 
ing. The model can be further simplified, however, 
by averaging over microion degrees of freedom to 
map the macroion-microion mixture onto a coarse- 
grained one-component model governed by effective 
electrostatic interactions [39] . The mapping acts on 
the partition function, 



Z = «exp(-/3iO) M ) 



(1) 



where H is the total Hamiltonian, f3 = l/k B T, and 
angular brackets denote traces over microion (/i) and 
macroion (m) degrees of freedom. The Hamiltonian 
naturally decomposes, according to H = Hm + H^ + 
H mll , into a bare macroion Hamiltonian H mi a mi- 
croion Hamiltonian H^, and a macroion-microion 
interaction energy H mfi . For a chemically closed 
suspension, which exchanges no particles with its 
surroundings, a canonical trace over only microion 
coordinates yields the canonical partition function 



Z = (cM-PH cS )) r 



(2) 



where H c g = H m +F^ is the effective one-component 
Hamiltonian and 

Ffi = -k B T\n (cxp + H mii )\) ^ (3) 

is the Helmholtz free energy of a microion gas in 
the midst of fixed macroions. Equations (2) and (3) 
provide a formal basis for approximating effective 
electrostatic interactions and simulating the effective 
one-component model of charged colloids. 



where -Fpiasma is the free energy of a uniform plasma 
of microions in a charge-neutralizing background 
of energy Eb, the charging parameter A tunes the 
macroion charge (and microion response) from zero 
to maximum, and ( ) A represents an average with 
respect to an ensemble of macroions charged to a 
fraction A of their full charge. For weakly correlated 
microions, the plasma free energy has the ideal-gas 
form, 

/plasma = iV+Hn+A 3 ) - 1] + 7V_[ln(n_A3) - 1], 

(5) 

where n± = N±/\V(1 — rj)\ are the average microion 
number densities, r\ = (47r/3)n m a 3 is the volume 
fraction of the macroions with number density 7i m — 
N m /V, and A M is the microion thermal wavelength. 

The linear-response approximation expands the 
microion number densities in functional Taylor se- 
ries in powers of the macroion external potential, 
truncates the expansions at linear order, and ne- 
glects microion correlations by assuming mean-held 
response functions [19-22]. The resulting internal 
potential energy, 

U = E vol (N m ,N s , V) + C/ pair ({r}; N m , N s , V), (6) 

separates into a one-body volume energy E vo \, 
which is independent of macroion coordinates, and 
a pair potential energy U pa i r , which depends on the 
macroion coordinates {r}. The volume energy, orig- 
inating from the microion entropy and macroion- 
microion interaction energy, is given by 



l3E vol = (3F pl 

a 



v,„ [ - 

z 



X, 



2 1 + KO 



Z n + — n_ 



(7) 
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where k = y/ 'A / n\BZ 2 n ll is the Debye screening con- 
stant (inverse screening length), a function of the 
total microion density, n M = n + + n_. The pair 
potential energy, 



t^pair — 2 



N m 

£ 



^off(|r 2 



(8) 



is a sum of hard-sphere-repulsive- Yukawa (screened- 
Coulomb) effective pair potentials, 



( cxp(/cq) \ 2 exp(-Kr) 
«eff(r) = <( yTT^a" J r > r ± 2a ' 

oo, r < 2a. 

(9) 

The effective pair potential, a product of microion 
screening of the bare macroion-macroion Coulomb 
interactions, has the long-range form of the DLVO 
potential [10, 11], but with a density-dependent 
screening constant. The constraint of electroneu- 
trality ties average macroion and microion number 
densities via Zn m /(1 — rj) = z(n + — n_), render- 
ing the effective interactions dependent on the aver- 
age densities of both macroions and salt ion pairs, 
n s = N S /[V(1 — rj)]. Equations (5)-(9) summarize 
the effective interactions that we input to theory and 
simulations of the one-component model. 



excess free energy density as 

/ex(n m ,n s ) = min \ fns{n m ,n s ; d) + 2^n 2 m 
(<*) I 

x J drr 2 g HS (r, n m ; d)v cS (r,n m ,n s )^ , (11) 

where the effective HS diameter d is the variational 
parameter and fns{ n m, n s] d) and <7Hs( r i n m\ d) are, 
respectively, the excess free energy density and (ra- 
dial) pair distribution function of the HS fluid, com- 
puted here from the near-exact Carnahan-Starling 
and Verlet-Weis expressions [38]. According to the 
Gibbs-Bogoliubov inequality [38], minimization of 
/ cx with respect to d yields a least upper bound to 
the free energy. From the variational approxima- 
tion for the total free energy [Eqs. (10) and (11)], 
the fluid branch of the phase diagram can be com- 
puted by performing a common-tangent construc- 
tion on the curve of free energy density f = F/V 
vs. macroion number density n m at fixed salt chem- 
ical potential, imposing equality of the pressure, 
P = n m(df /dn m ) Ng / Nm — /, and of the macroion 
and salt chemical potentials, \x m = (df /dn m ) ns and 
ji s = (df /dn s ) nm , in coexisting phases. 



IV. MONTE CARLO SIMULATIONS 



B. Variational Free Energy Theory 

At constant particle numbers, volume, and tem- 
perature, the Hclmholtz free energy F is a minimum 
with respect to variations in N m , N±, V, and T 
at thermodynamic equilibrium. The electroneutral- 
ity constraint requires that ion exchange between 
phases occurs only in charge-neutral units, allowing 
the free energy to be regarded as a function of the 
number of salt ion pairs N s , rather than of N + and 
iV_ separately. Within the one-component model, 
the free energy separates, according to 

F(N m ,N s ,V) - F id (N m ,V) + F ex (N m ,N s ,V) 
+ E vol (N m ,N s ,V), (10) 

where F^ = iV m fcsT[ln(n m A 3 ) — 1] is the free en- 
ergy of an ideal (noninteracting) gas of macroions 
of thermal wavelength A, and F ox is the excess free 
energy due to effective pair interactions [Eq. (9)]. 

A variational approximation [12-14, 23] based on 
first-order thermodynamic perturbation theory with 
a hard-sphere (HS) reference system [38] gives the 



The effective interactions described above, which 
were used in previous variational theory calculations 
for the one-component model [23], are here input 
into simulations of the same model to test the ac- 
curacy of the variational approximation and its pre- 
dictions for thermodynamic behavior. The Gibbs 
ensemble Monte Carlo (GEMC) method [40-44] is 
an efficient means of simulating two-phase fluid co- 
existence that obviates the need to model interfaces. 
Each phase is represented by its own simulation 
box, with fluctuating macroion numbers N m i and 
volumes Vi (i — 1,2). In the constant-AVT im- 
plementation, the total macroion number, N m = 
N m \ + N m 2, total volume, V = V\ + V2, and tem- 
perature T all remain fixed. We further fix the total 
number of (implicit) salt ion pairs, N s = N s \ + N S 2, 
while performing virtual exchanges between boxes. 
Although the GEMC method has been previously 
applied to fluids with density-dependent pair poten- 
tials [45], it has not yet, to our knowledge, been 
adapted to charged systems whose effective interac- 
tions include both a pair potential and volume en- 
ergy. 

The conventional GEMC algorithm [40-44] in- 
volves three types of random trial move: (1) dis- 
placements of particles (macroions) within each box 
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to ensure thermal equilibrium of each phase; (2) 
volume exchanges between the two boxes to ensure 
mechanical equilibrium, characterized by equality of 
pressures; and (3) macroion transfers between the 
two boxes to ensure chemical equilibrium with re- 
spect to macroion exchange, characterized by equal- 
ity of macroion chemical potentials. The acceptance 
probability P movc for any trial move from an old (o) 
to a new (n) state can be derived from the Metropo- 
lis condition [46-48], 



P move = min <! t^4, 1 



V(o) 



(12) 



where the Gibbs ensemble probability density [47] is 
given by 



V = 



1 



N ml \N m2 \ v as ; v ,v> 

x exp[-{3U({s};n m ,n s )} 



(13) 



and {s} denotes the macroion coordinates scaled by 
their respective box lengths. Although the salt ion 
coordinates do not explicitly appear in Eq. (13), the 
potential energy U [Eq. (6)] implicitly depends on 
the average salt (and macroion) densities in the two 
boxes. 

From Eqs. (12) and (13), trial displacements are 
accepted with probability 



Pdisp = min {cxp(-/3AP), 1} . 



(14) 



where AU = U(n) — U(o) is the change in total 
potential energy between the new and old states. 
Note that for internal displacements, which do not 
affect the volume energy, AU = A[/ pa i r [Eq. (8)]. 
For all other moves, however, the change in total 
potential energy also includes a change in volume 
energy: AU = AE vol + AU pail . 

A trial exchange of volume AV from box 1 to 
box 2 (Vi -> Vi - AV, V 2 ^V 2 + AV) is achieved 
by uniformly rescaling all macroion coordinates. In 
practice, it proves more efficient to vary ln(Vi/V2), 
with an acceptance probability [47] 



P vo i = min' 



■V-AV^^ 1 /V 2 + AV^ Nm2+1 



v 2 



x exp(-/3A£/), 1 \. 



(15) 



Transfer of a macroion from box 1 to box 2 {N m \ — > 
N m i ~ 1, N m2 — > N m2 + 1) is accepted with proba- 



bility [42] 



P 



trans 



N„ 



V 



N m2 + 1 V 1 



exp(-/3AP), 1 



(16) 

Note that AU in Eqs. (15) and (16), represents the 
change in total potential energy of the two boxes 
combined, since exchanges of volume or macroions 
alter the average macroion density (n m i — N m i/Vi), 
and thus the volume energy [Eq. (7)] and pair po- 
tential [Eq. (9)], in each box. 

In addition to the conventional GEMC moves, we 
introduce a new trial move: transfer of salt between 
the two boxes, required to ensure chemical equilib- 
rium with respect to salt exchange between coexist- 
ing phases, characterized by equality of salt chemi- 
cal potentials. Since the salt is modeled here only 
implicitly, virtual transfers involve simply changing 
the average salt density of each box, with acceptance 
probability 



P sa it =min{cxp(-/3AP),l}. 



(17) 



where AU is the change in total potential energy 
of both boxes. We stress that exchanges of average 
salt density affect both the pair potential and the 
volume energy in each box. The absence of combi- 
natorial and phase-space prefactors in Eq. (17) fol- 
lows from implicit modeling of salt ions. In practice, 
a transfer of AiV^ salt ion pairs from box 1 to box 2 
(iV s i N sl - AN S , N s2 -> N s2 + AN S ) is realized 
by changing the respective salt densities accordingly 
and adjusting AN S N s i,N s2 ) to achieve a rea- 
sonable acceptance rate. 

Within the Gibbs ensemble, we simulated two cu- 
bic boxes subject to periodic boundary conditions, 
each box containing only macroions, but evolving 
according to effective interactions [Eqs. (7) and (9)] 
that implicitly depend on the microion densities. To 
exclude interactions of a particle with its own peri- 
odic images, and avoid needless computation, pair 
interactions between macroions were cut off at a dis- 
tance of r c = min{20/K, L/2}, i.e., the shorter of 20 
screening lengths or half the respective box length 
L. The effective interactions were updated whenever 
the average macroion or salt density changed. 

The simulations started from initial configurations 
of randomly distributed macroions, with equal par- 
ticle numbers, volumes, and salt concentrations in 
each box. The four types of trial move were ex- 
ecuted in random sequence at prescribed frequen- 
cies. Defining a cycle as an average of N m trial dis- 
placements (i.e., one per macroion), the other moves 
were attempted with relative frequencies per cycle 
of N m /2 for volume exchanges, N m /10 for macroion 
transfers, and N m /10 for salt exchanges. For inter- 
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nal displacements, macroions were selected at ran- 
dom and moved with tolerances adjusted to yield 
an acceptance rate of about 50%. For volume and 
salt exchanges, acceptance rates of about 10% were 
achieved by adjusting the tolerances, the resulting 
salt tolerance being AN S /N S ~ 10~ 3 . After equili- 
brating for 10 4 cycles, we accumulated statistics for 
average densities, pressures, and chemical potentials 
over the next 10 4 cycles (5 x 10 6 displacements for 
N m = 500). 



V. RESULTS AND DISCUSSION 

A. Tests of One-Component Model 
and Variational Theory 

To investigate thermodynamic phase behavior of 
charged colloids, we input effective electrostatic in- 
teractions (Sec. Ill A) to both variational theory cal- 
culations (Sec. IIIB) and Gibbs ensemble Monte 
Carlo simulations (Sec. IV) of the coarse-grained 
one-component model. The validity of the one- 
component model is first tested by comparing struc- 
tural and thermodynamic properties with available 
data from simulations of the primitive model, which 
include explicit point counterions interacting via 
bare Coulomb potentials. 

From extensive Monte Carlo (MC) simulations, 
Linse [35] has generated a wealth of data for the 
(salt-free) primitive model over ranges of macroion 
valence, volume fraction, and electrostatic cou- 
pling parameter, T = As/a. For direct compari- 
son, we performed simulations of the effective one- 
component model for identical parameters - fixing 
the effective macroion valence (Z = 40), counterion 
valence (z = 1), and Bjerrum length (As = 0.72 
nm), and varying the macroion radius a - and com- 
puted the macroion-macroion pair distribution func- 
tion g(r) and pressure p } as described in the Ap- 
pendix. For this purpose, we performed standard 
constant-./WT (one-box) simulations, the volume 
energy then having no effect on the pair structure. 
To obtain accurate g(r) results, a system size of 
N m = 600 sufficed to render finite-size effects negli- 
gible. To maintain consistently high accuracy in the 
pressure, we increased the particle number to en- 
sure a cut-off radius of at least 20 screening lengths 
for each combination of -q and T — ranging up to 
N m = 8000 for r, = 0.02, T = 0.0445. 

Figures 1 and 2 compare numerical results of 
our simulations of the one-component model with 
Linse's simulations of the primitive model [35]. At 
relatively low electrostatic couplings (r < 0.1779, 
ZXb/cl < 7), our results for the pair distribution 



function and pressure closely match the correspond- 
ing primitive model data (Fig. 5 (a) and Table III 
of rcf. [35], after minor corrections [49]). It should 
be noted that good agreement at higher volume frac- 
tions is achieved only when the excluded- volume fac- 
tor of 1/(1 — 77) is consistently included in the effec- 
tive interactions. These comparisons demonstrate 
the accuracy of the one-component model with lin- 
earized effective interactions for moderately coupled 
systems. Figure 2 also presents predictions for the 
pressure from our variational theory calculations. 
The near-perfect alignment of theory and simula- 
tions of the one-component model validates the vari- 
ational approximation over the parameter ranges 
studied. 
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FIG. 1: Macroion-macroion pair distribution function 
g(r) vs. radial distance r (units of macroion radius a) of 
salt-free suspensions computed from Monte Carlo sim- 
ulations of the one-component model (OCM) with im- 
plicit counterions (solid curves) and the primitive model 
(PM) [35] with explicit counterions (dashed curves) for 
macroion valence Z = 40, volume fraction r\ = 0.01, and 
electrostatic coupling parameters F = As/a = 0.0222, 
0.0445, 0.0889, 0.1779, 0.3558 (bottom to top). For clar- 
ity, curves are vertically offset in steps of 0.5. 

At higher electrostatic couplings {ZXb/cl > 7), 
typical of highly charged latex particles and ionic 
surfactant micelles, significant deviations between 
the one-component and primitive models abruptly 
emerge (r = 0.3558 in Figs. 1 and 2). The discrepan- 
cies in this relatively strong-coupling regime can be 
traced to renormalization of the effective macroion 
charge through strong association of counterions, 
a nonlinear effect neglected in the present version 
of the model. Preliminary investigations [50] indi- 
cate, however, that the deviations can be substan- 
tially reduced by consistently building into the one- 
component model a renormalized effective charge. 
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FIG. 2: Reduced pressure /3p/n to t vs. macroion volume 
fraction r/, where n to t = (Z + l)n m , for salt-free sus- 
pensions computed from Monte Carlo simulations of the 
effective one-component model with implicit counterions 
(solid symbols) and the primitive model [35] with explicit 
counterions (open symbols) for macroion valence Z = 40 
and several electrostatic coupling parameters T = As /a. 
Simulation error bars are smaller than symbol sizes. Also 
shown are corresponding predictions of variational the- 
ory (curves). From top to bottom, T = 0.0222, 0.0445, 
0.0889, 0.1779, 0.3558. 



These results establish a threshold of Z\s/a ~ 7 
for significant charge renormalization within linear- 
response theory. 

To test the variational free energy theory at higher 
charge asymmetries and nonzero salt concentrations, 
we compare predictions for the osmotic pressure 
(equation of state) with results from our GEMC sim- 
ulations. The osmotic pressure, II = p — 2n r kBT, is 
here defined as the total pressure of the suspension 
less that of a (virtual) ideal-gas salt reservoir of ion 
pair density n r at the same salt chemical potential: 
l^s = kBT\n(2n r A^). Figure 3 shows sample results 
for the equation of state at fixed salt chemical poten- 
tial or, equivalently, salt fugacity, z s — exp(/3/z s /2), 
from simulations at the salt concentrations predicted 
by theory for each volume fraction. Since theory and 
simulation assume identical effective interactions, 
the comparisons directly probe the excess free energy 
approximation [Eq. (11)] and corresponding pair po- 
tential contribution to the total pressure (inset to 
Fig. 3). The predictions are in excellent agreement 
with simulation over a wide range of volume frac- 
tions, further validating the variational approxima- 
tion and providing a consistency check on our calcu- 
lations. As an independent check, our methods accu- 
rately reproduce pressures computed from MC simu- 
lations of the hard-sphere-repulsive- Yukawa pair po- 
tential fluid [51]. 

The appearance in Fig. 3 of a van der Waals loop 
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Macroion Volume Fraction r\ 

FIG. 3: Reduced osmotic pressure /3IIa 3 vs. volume frac- 
tion r] computed from Monte Carlo simulations and vari- 
ational theory [23] of the one-component model for Bjer- 
rum length As = 0.72 nm, macroion radius a = 50 nm, 
valence Z = 500, and salt fugacities z s — 100 and 200 
/jM. Changes of curvature reflect phase instability. Inset: 
Pair potential contribution to total pressure. Simulation 
error bars are smaller than symbol sizes. 

in the pressure signals a spinodal instability and sep- 
aration into macroion-rich (liquid) and macroion- 
poor (vapor) phases. We stress, however, that cur- 
rently available data from primitive model simula- 
tions can test the effective one-component model 
and linearized effective interactions only for salt- 
free suspensions at relatively low charge asymme- 
tries, where instabilities with respect to phase sepa- 
ration have not been predicted. Furthermore, the 
macroion aggregation observed in ref. [35] in the 
strong-coupling regime is likely driven by microion 
correlations, which are neglected in the mean-field 
effective interactions assumed here. While further 
tests of the one-component model are needed, the 
close agreement for parameters accessible to primi- 
tive model simulations motivates proceeding to con- 
sider phase behavior. 



B. Phase Behavior 

To systematically map out the fluid binodal, we 
performed a series of GEMC simulations over ranges 
of volume fraction and salt concentration for selected 
macroion radii and valences: (a — 10 nm, Z = 150) 
and (a = 50 nm, Z = 500). Initially uniform systems 
of N m = 500 particles (two-box total), in thermody- 
namic states (77, c s ) within the predicted unstable re- 
gion [23], spontaneously separated into two phases, 
each phase occupying one of the boxes, which dif- 
fered in average macroion and salt densities. In con- 
trast, systems at state points outside of the unstable 
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region remained uniform. A visual impression of the 
phase separation is provided by the simulation snap- 
shot in Fig. 4. 



quantitative agreement between theory and simula- 
tion attests to the accuracy of the variational ap- 
proximation. 




FIG. 4: Snapshot from Gibbs ensemble Monte Carlo 
simulation, showing the two boxes after separation into 
macroion-rich and macroion-poor phases. Spheres depict 
macroions in the effective one-component model. 

To identify the structure of the coexisting phases, 
we performed constant- NVT (one-box) simulations, 
at identical state points, for particle numbers com- 
mensurate with likely crystal structures: fee (N m — 
500) and bee (N m = 432). Initializing the particles 
on the sites of the respective lattice, we computed 
the equilibrium pair distribution function and ob- 
served typical fluid-like structure, indicating melt- 
ing of the initial crystal. Upon increasing the vol- 
ume fraction, we observed, at state points well out- 
side the fluid binodal, an abrupt sharpening of the 
peaks of g(r), reflecting crystallization. These ob- 
servations are consistent with a simple hard-sphere 
freezing criterion, rj(d/2a) 3 ~ 0.49, which approxi- 
mates the macroions as hard spheres of effective di- 
ameter d [from Eq. (11)] and locates the coexistence 
densities within the fluid regime. 

The resulting phase diagrams are presented in 
Fig. 5, alongside predictions of variational the- 
ory [23], where tie lines joining corresponding points 
on the macroion-rich and macroion-poor binodal 
branches parallel those predicted by theory. Each 
pair of points on the binodal was produced by aver- 
aging over 10 independent runs, which differed only 
in the random number seed used for trial moves. Re- 
ported error bars represent statistical uncertainties 
of one standard deviation, computed from fluctua- 
tions in average densities among the 10 runs. Res- 
olution near the critical point is blurred by density 
fluctuations and phase switching between boxes - a 
known limitation of the Gibbs ensemble method [47] . 
For simplicity, we discarded runs in which the phases 
switched boxes, a rare occurrence away from the crit- 
ical point. Considering the sensitivity of the coexis- 
tence analysis to slight deviations in free energy, the 




400 



o 300 



c 
o 

"-i — • 

CO 
k_ 
-t— ' 

c 

CD 

o 
c 
o 
O 

-i — ' 

CO 




0.01 0.02 0.03 0.04 0.05 0.06 

Macroion Volume Fraction r\ 

FIG. 5: Phase diagrams of aqueous charged colloids, 
showing fluid binodal computed from Gibbs ensemble 
Monte Carlo simulations (squares) and variational the- 
ory [23] (solid curves) for the effective one-component 
model at various combinations of macroion radius a and 
valence Z. Tie lines join corresponding points on the 
colloid-rich and colloid-poor branches of the binodal. 
Also shown are predictions of variational theory for the 
critical point (circles) and spinodal (dashed curves). 

Diagnostic variables were monitored during the 
simulations and evolved as typified by Fig. 6, which 
tracks the volume fractions, salt concentrations, 
pressures, and chemical potentials in each box vs. 
number of MC cycles for one sample run. Bifur- 
cation of volume fractions and salt concentrations 
in the two boxes [Fig. 6 (a)] signals phase separa- 
tion, while convergence of pressures and chemical 
potentials [Fig. 6 (b) and (c)] confirms equilibra- 
tion. Several runs for larger systems (up to 1000 
macroions) were performed to establish the insignif- 
icance of finite-size effects. Compared with large- 
scale simulations of the primitive model, our Simula- 



tions require relatively modest computing resources, 
each run typically consuming 50-90 CPU hours on a 
single PC (Intel Xeon-HT processor). 

It should be emphasized that the observed phase 
separation, although perhaps surprising in the face 
of purely repulsive pair interactions, is driven by the 
state dependence of the volume energy in the one- 
component model of deionized suspensions. It is also 
important to note that the macroion valences and 
electrostatic couplings represented in Fig. 5 were se- 
lected to lie just within the unstable fluid regime 
and correspond to T = 0.072 and 0.014 (ZX B /a = 
7 and 11) in panels (a) and (b), respectively (c/. 
ZX B /a ~ 1-14 in ref. [35]). In each case, a small 
increase in macroion radius or decrease in valence 
stabilizes the system. These parameters approach 
the threshold for charge renormalization estimated 
from our direct comparisons with primitive model 
simulations (Figs. 1 and 2), albeit for lower valences. 
Whether the predicted phase instability corresponds 
to a real phenomenon or is merely an artifact of 
the linearization approximation [27-29], or the as- 
sumption of fixed macroion valence [30-34] , remains 
unclear. Preliminary explorations [50], based on a 
simple model of effective macroion valence, suggest 
that the instability survives incorporation of charge 
renormalization in the one-component model. Fur- 
ther studies are required, however, to resolve this 
important issue. 



VI. CONCLUSIONS 

In summary, we have developed a new variant of 
the Gibbs ensemble Monte Carlo method to simulate 
a one-component model of charged colloids governed 
by density-dependent effective interactions. The ef- 
fective interactions (pair potential and one-body vol- 
ume energy) are input from linear-response theory, 
assuming a mean-field approximation for microion 
structure. The simulation algorithm includes trial 
exchanges of implicit salt between the two simu- 
lation boxes and incorporates the volume energy 
into the acceptance probabilities for trial moves that 
change the average macroion or salt density. 

Comparisons with simulations of the two- 
component (salt-free) primitive model [35] demon- 
strate the validity of the one-component model 
over a wide parameter range, physically relevant to 
charged latex particles and micelles. Results for the 
macroion-macroion pair distribution function and 
pressure are in close agreement with corresponding 
primitive model results for moderate electrostatic 
couplings. Deviations at stronger couplings likely 
originate from nonlinear screening effects neglected 
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FIG. 6: Evolution of diagnostic properties during a sam- 
ple Gibbs ensemble Monte Carlo run for macroion ra- 
dius a = 50 nm, valence Z = 500, total volume fraction 
r\ — 0.025, and total salt concentration c s — 200 /iM. 
Solid and dashed curves represent the two boxes, (a) vol- 
ume fractions and (inset) salt concentrations; (b) total 
pressure; (c) macroion and (inset) salt chemical poten- 
tials, shifted to zero average. Instantaneous values are 
plotted during equilibration (first 10 4 cycles) and cumu- 
lative values thereafter. 



in the present model. Our simulations also confirm 
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the accuracy of a variational free energy approxima- 
tion [12-14, 23]. Further comparisons would help 
to more sharply define the limitations of the one- 
component model. 

While the cost of primitive model simulations 
grows with increasing charge asymmetry, concentra- 
tion, and electrostatic coupling, the computational 
effort required to simulate the one-component model 
is relatively modest and actually diminishes with 
increasing macroion valence and concentration, as 
decreasing the screening length shortens the range 
of effective pair interactions. The one-component 
model thus offers insight into bulk phase behavior 
in parameter regimes that may be computationally 
prohibitive for more explicit models. 

We have applied our new simulation method to 
test predictions of variational theory [23] for the 
phase behavior of aqueous suspensions of charged 
macroions with weakly correlated (monovalent) mi- 
croions at low salt concentrations. The resulting 
phase diagrams exhibit coexistence of macroion-rich 
and macroion-poor fluid phases in excellent agree- 
ment with previous predictions and qualitatively 
consistent with observed thermodynamic anomalies. 
The phase instability predicted by theory, and now 
confirmed by simulations of the same model, occurs 
in a parameter regime that appears to border the 
threshold for saturation of the effective macroion 
charge. Future work will address this open issue 
by incorporating charge renormalization into the 
one-component model [50]. Finally, our simula- 
tion algorithm can be extended to investigate other 
phase transitions, e.g., crystallization, and adapted 
to model other soft materials, such as polyelectrolyte 
and ionic micellar solutions. 
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APPENDIX: DIAGNOSTIC CALCULATIONS 

1. Pressure 

A diagnostic for mechanical equilibrium in the 
Gibbs ensemble is equality of pressures in the two 
boxes. The total pressure naturally separates into 
three distinct contributions: 



where p;d — (n m ) ksT is the ideal-gas pres- 
sure of the macroions, p pa j r results from effec- 
tive pair interactions between macroions, p vo \ = 
— ((dE vo \/dV)]y m /N s ) is the contribution from 
the density-dependent volume energy, and angular 
brackets denote an ensemble average over configu- 
rations in the Gibbs ensemble. The pair pressure is 
calculated on the fly within the simulations using the 
virial expression for a density-dependent pair poten- 
tial [52]: 



Ppa 




l+ptail, (A.2) 



N m /N„ 



where V; n t is the internal virial, the volume deriva- 
tive term accounts for the density dependence of the 
effective pair potential, and p ta ii corrects for cutting 
off the long-range tail of the pair potential. The 
internal virial is given by 

N m N m 
Vint = ^ *i ■ fj = ^2 ( 1 + Kr ij) V cS(rij), (A.3) 

i=l i^j=l 

where fj = — V J2j^i v cS ( r ij) 1S the effective force 
exerted on macroion i, at position r i; by all neigh- 
boring macroions j, at relative distances , within 
the cut-off radius r c . The second term on the right 
side of Eq. (A.2) is computed via 



dV ) N s /N m ~ 2Nm t ^tl ^ ;hi ' ! 



E 



with 



N 3 /N m 

(A.4) 



/ dv cS (r) ' 

V Q n m / N s /N„ 



k a nr\ v e g(r) 



1 + na 2 J n m (l - r]) ' 
(A.5) 

The tail pressure is approximated by 



Ptaii = _ y( n m^ drr 3 v' eS (r) 



2tt 



K 2 r 2 r + 3nr c + 3 



r c v e s(r c ) ) , 



(A.6) 



the approximation being the neglect of pair corre- 
lations for r > r c . Finally, the volume pressure is 
given by 



1-T) 



Z + 2 



_ Z 2 k\ b 
N m 4(1 + na) 2 



(A.7) 



P = Pid +Ppair +Pvol, 



(A.l) 
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2. Chemical Potentials 

To diagnose chemical equilibrium between coexist- 
ing phases, we computed the chemical potentials of 
macroions and salt by adapting Widom's test parti- 
cle insertion method [53] to the Gibbs ensemble, fol- 
lowing ref. [47]. In contrast to the original method, 
the inserted ions are not treated as ghost particles in 
GEMC, but rather remain within the box into which 
they are successfully transferred. The macroion 
chemical potential - the change in Helmholtz free 
energy upon adding a macroion - can be expressed 
as 

^ V Z G (N m ,N s ,V,T) )' y ' 

where the Gibbs ensemble partition function is given 

by 



z G 



l 



N„ 

E 



i 



A3NmV ivflto N "* l ( N ™ ~N ml )l 



^N m -N„ 



x f dVtV^iV-Vt) 1 
Jo 

x Jds N ™ exp[-/?[/({s};n ro ,n s )]. (A.9) 

The macroion chemical potential of box 1 is thus 
computed from [47] 



Mmi = -k B T In 



V 1 /A 3 



■exp(-/?At/+™)), 



, (N ml + 1) 

(A.10) 

where AUf™ = U(N ml + 1) - U(N ml ) is the change 
in total potential energy (volume energy plus pair 
energy) of box 1 upon insertion of a macroion. In 
practice, the large change in volume energy AE vo i 
resulting from a macroion insertion necessitates eval- 
uating Eq. (A. 10) by adding to and subtracting 
from the argument of the exponential a constant 



Ihni = -k B Tln( ^ l/A " exp[-/?(At/+' 

+ C. 



(ATI) 



The salt chemical potential - the change in 
Helmholtz free energy upon insertion of a salt ion 
pair - can be approximated by 

k B T ( Z G (N m ,N s + AN s ,V,T) 
^ AN S U { Z G (N mi N s ,V,T) 

(A.12) 

assuming that the number of exchanged salt ion 



pairs AN S is much less than the total number of salt 
ion pairs (AN S <C N s ). The salt chemical potential 
of box 1 is thus computed from 



^ 1 = ~A^ ln ^ exph/3(AC/l+S ~ c) ^ 



AAV 
(A.13) 

where AU^ S = U(N sl +AN s )-U(N sl ) is the change 
in total potential energy of box 1 upon insertion 
of AN S salt ion pairs. The absence of combinato- 
rial and phase space factors in Eq. (A.13) follows 
from modeling the microions only implicitly. Note 
also that the chemical potentials are defined only 
to within arbitrary constants. In the dilute colloid 
limit (N m — > 0), the salt chemical potential tends to 
that of an ideal gas of salt ions 



A 4 . 



(o) 



2fc B T<ln(n s A3)> 



(A.14) 



and the macroion chemical potential reduces to 



(o) _ 



= k B T I - In 



Z\n{n s Kl) 



: / 2(1 + Kfl ) + T^ 3 )' (A - 15) 



where the terms on the right side are derived (left to 
right) from the macroion entropy, microion entropy, 
macroion-counterion interaction, and macroion ex- 
cluded volume. These analytical results [Eqs. (A.14) 
and (A. 15)] provide a check on the numerical results 
in the limit in which one box becomes depleted of 
macroions. 



3. Pair Distribution Function 

The structure of the suspension is character- 
ized by the pair distribution functions [38]. The 
macroion-macroion pair distribution function g(r) 
- the only one accessible in the one-component 
model - is defined such that Airr 2 g(r)dr equals 
the average number of macroions in a spherical 
shell of radius r and thickness dr centered on a 
macroion. For a given configuration, each particle 
is regarded, in turn, as the central particle. Neigh- 
boring particles are then assigned, according to 
their radial distance r from the central particle, 
to concentric spherical shells (bins) of thickness 
Ar = 0.1a. After equilibration, g(r) is computed, 
in the range 2a < r < L/2, by accumulating the 
numbers of particles in radial bins and averaging 
over all configurations. The resulting distributions 
are finally smoothed by averaging each bin with its 
immediate neighboring bins. 
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